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Abstract. When a magnetised, thermal plasma containing a 
shock wave efficiently accelerates Cosmic Rays (CR) a reac- 
tion is exerted on the fluid structure. The simplest description 
of this process involves three conservation equations for the 
^\ background fluid, with the CR pressure gradient included as 
^\ an extra force, plus a hydrodynamic equation for the CR en- 
ergy density. Since the CR spectrum is not determined in this 
^jf^approach two quantities are introduced, the ratio of CR pres- 
1-^ sure to energy density and a mean diffusion coefficient, that 
^rfH must be specified to close the system of equations. Physical ar- 
guments are presented for how these quantities should evolve 
and then compared with time dependent numerical solutions 
of acceleration at a plane parallel, piston driven shock. 

T— H Key words: acceleration mechanisms - cosmic rays - shock 
^ 1 ^ waves - magnetohydrodynamics - plasmas - supernova rem- 
00 nants 
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Introduction 

O ^There are a variety of astrophysical environments where the 
I energy density of Cosmic Rays (CRs) is comparable to that of 
O the local thermal plasma. Such is the case, for example, in the 
interstellar medium (ISM) where the dynamics of the plasma is 
^ coupled to that of the energetic particles not only through the 



latter's pressure and energy flux but also by collective wave- 
particle excitation. The combined dynamics of thermal and 
non-thermal particles is important not just for the propaga- 
tion of CRs in the ISM but also during the CR acceleration 
process itself which is believed to occur, at least up to the knee 
of the CR spectrum, at the strong shock waves of supernova 
remnants (SNRs) by the first order Fermi mechanism (Axford 
et al. 1977, Krymsky 1977, Bell 1978 a,b and Blandford and 
Ostriker 1978). The energetic particles are scattered by the 
magnetic irregularities which are essentially frozen into the up- 
stream and downstream fiows and gain energy by compression 
at the shock front. Once the pressure of the accelerated CRs 
rises to a significant fraction of the gas pressure their effect on 
the hydrodynamics becomes important so that a pure thermal 
gas dynamic approach to the evolution of the SNR is no longer 
accurate. The simplest way of calculating the hydrodynamic 
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modification is to include the CR pressure, Pc, as an extra 
component in the total pressure. Staying with the example of 
a SNR, with the region ahead of the expanding shock front 
occupied by a diffusive tail of CRs the upstream plasma is ac- 
celerated, by this CR pressure gradient, and possibly heated, 
through the dissipation of resonantly excited Alfven waves. The 
expanding blast wave then sweeps up not the undisturbed ISM 
nor the progenitor star's stellar wind but a mixture of a modi- 
fied thermal plasma and non-thermal CRs. Furthermore since 
the energetic particles are scattered by the magnetic turbu- 
lence there is a diffusive coupling between the modified up- 
stream precursor and the high pressure, low density interior of 
the remnant. With efficient particle acceleration the explosion 
energy of a supernova is therefore converted in a non-linear 
fashion not just into kinetic and thermal gas energies but also 
into non-thermal CR energy (Drury et al. 1989, Dorfi 1990, 
Jones and Kang 1990, Markiewicz et al. 1990 and Berezhko et 
al. 1993). This model of self consistent SNR and CR evolution 
can provide, with reasonable choices of parameters, both suffi- 
cient power to sustain the observed CR intensity and enough 
hot gas to explain the X-ray luminosity of young SNRs. 

Crucial to the above picture, however, is exactly how the 
CR energy density, Ec, and pressure are to be calculated. In 
general these bulk quantities depend on an energetic particle 
spectrum that is neither in equilibrium nor a simple power law, 
since the shock structure is modified. This suggests that one 
should solve self consistently for the phase space density dis- 
tribution function f[x,p,t) of the energetic particles and from 
this obtain the fiuid moments Ec and Pc (Falle and Giddings 
1987, Bell 1987, Duffy 1992 and Berezhko et al. 1993). An 
alternative, and numerically simpler, approach which is sug- 
gested intuitively by the notion of a bulk CR pressure and 
energy density is to extend the gas dynamic view to include 
the CR population without obtaining information about the 
spectrum (Drury and Volk 1981 and Axford et al. 1982). The 
almost isotropic distribution of CRs is treated as a fiuid with 
a negligible mass density but an important, and in some cases 
dominant, energy density. In contrast to the thermal energy 
density which, in a comoving volume element with dissipative 
heating ignored, can only change as a result of adiabatic com- 
pression the added effect of scattering off magnetic turbulence 
is included for Ec- While this picture of a compressible, dif- 
fusing fiuid which provides an extra component to the total 
pressure is straightforward qualitatively, agreement with a ki- 
netic solution of the CR transport is not so simple. Turning 
to compression first, as the volume element is deformed the 
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Fig. la Modified shock profile with Uo, Ui and U2 the far up- 
stream, pre-shock and downstream flow speeds respectively. 

change in Ec will depend on the CR pressure given by the 
equation of state, Pc = {fc — i)Ec where fc is the CR 'adi- 
abatic index'. During the acceleration process fc will evolve 
from 5/3 to 4/3 as ultrarelativistic particles begin to dominate 
the spectrum. Secondly, since on a kinetic level the CR diffu- 
sion will be momentum dependent with a coefficient /c(p), the 
rate of diffusion of the bulk energy density will be given by 
some mean 7c of the microscopic coefficient weighted over the 
CR spectrum (see equation 3.2 below). Therefore CR hydro- 
dynamics at shock fronts has the difficulty that the parameters 
fc and 7c depend on the evolving spectrum of energetic parti- 
cles which is undetermined; the closure problem. Since, for the 
case of acceleration at a shock front, the solutions to the two 
fluid model are quite strongly dependent on the choice made 
for fc and 7c (Achterberg et al. 1984 and Heavens 1984) the 
challenge is to find physical arguments that would allow us to 
predict the evolution of the closure parameters in terms of fiuid 
quantities. 

The first attempt in this direction was contained in Drury 
et al. (1989) in connection with the simplified shell models. 
The expanding supernova remnant is divided into a number 
of regions with conservation laws describing the dynamics of 
each. Further work by Kang and Drury (1992) showed that 
these simplified and two fiuid models agree reasonably well in 
their predictions. On the basis of sub-shock advection and in- 
jection Drury et al. estimate what proportion of the total CR 
energy density is made up of non-relativistic particles. Assum- 
ing that the remainder are ultra-relativistic this gives fc as 
a weighted mean of 5/3 and 4/3. The mean diffusion coeffi- 
cient is calculated by first obtaining the time dependent upper 
cut-off to the spectrum, Pmax. A fraction of /c(pmax) is then 
taken for 7c. Two points can be made about these arguments. 
Firstly the closure parameters depend on the CR spectrum so 
that it would be more direct to relate fc and 7c to anything 
that could be said about modified CR distributions. Second 
is that any such arguments should be tested against a full ki- 
netic code where available. It was not until Volk et al. (1991) 
(hereafter referred to as VDDM) that an attempt was made to 
construct a model with these points in mind. Use was made of 
the fact that in all but the completely cosmic ray dominated 
case, a sharp gas shock is modified to give a weaker subshock 



Fig. lb The double power law spectrum with Qs and qt given 
by the subshock and total compression ratios. 

followed by a precursor in the upstream region caused by the 
CR pressure gradient (figure la). With a diffusion coefficient 
increasing with energy the spatial scale of this precursor is de- 
termined by the diffusion lengthscale of a particle with some 
momentum, p., intermediate between the point of injection, pi, 
and the time dependent upper cut-off, pmax (eq. 8 of VDDM). 
Consequently while those particles near injection only feel the 
compression due to the gas subshock, the highest energy CR 
see the entire subshock-precursor structure as a discontinuity. 
In accordance with test particle theory it was assumed that 
the spectral slopes at pi and pmax are then determined by the 
subshock and total compression ratios respectively. With each 
of these slopes interpolated to the intermediate momentum p, 
this gives a double power law model for the modified CR dis- 
tribution function from which the closure parameters may be 
calculated (figure lb). Jones and Kang (1992) and Kang (1993) 
have used the time evolution of a test particle distribution func- 
tion to calculate the closure parameters using equation (38) of 
Drury et al. (1989). 

The double power law idea can be tested in the planar 
case through a direct comparison with a numerical kinetic so- 
lution to the problem. This is carried out below in section (2) 
where care is taken to consider a case for which the prediction 
of the closure parameters ought to be difficult. To this end 
we consider the limit of quasilinear theory, the Bohm limit, 
where the mean free path of a CR particle is equal to its gy- 
roradius in the mean magnetic field. This ought to be the case 
near strong shocks where wave self-excitation may be so strong 
that diffusion approaches this lower bound of /c <x pv, with v 
the particle's velocity. Therefore 7c will be a non-trivial av- 
erage of a strongly momentum dependent quantity weighted 
by a spectrum extending over several decades of energy. Sec- 
ondly we consider shock waves that are substantially modified, 
though not smoothed out completely, by the non-thermal pop- 
ulation. Finally in a planar shock with losses ignored, accel- 
eration will continue indefinitely so that pmax enters the ul- 
trarelativistic regime and fc will come close to 4/3. We there- 
fore restrict ourselves to intermediate timescales, i.e. timescales 
where most of the CR energy density is due to particles in the 
trans-relativistic regime so that fc is not close to its lower or 
upper limit. It will be found (section 2) that for strongly mod- 
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ified shocks in the Bohm limit during intermediate timescales 
that although the double power law approximation holds for 
the lowest energy CRs and although the spectrum flattens rel- 
ative to the subshock compression at high energies the discrep- 
ancy between the model and the computed spectra is too great 
for the former to accurately predict the closure parameters. 
However, there is recent evidence (Berezhko et al. 1993 and 
E.Berezhko, private communication) that for long timescales, 
beyond the timescales considered here, the CR spectrum does 
relax to a double power law and that the discrepancies we find 
below are endemic to our relatively small dynamical range. Dis- 
cussion of these results is left to the end of this paper after we 
have addressed what appears to be the most difficult part of the 
closure problem; accurately predicting fc and 7c for interme- 
diate timescales (section 3). Comparison with the steady state 
Monte Carlo simulations due to Ellison and Eichler (1984) is 
less fruitful because of the intrinsically time dependent char- 
acteristics of the problem under consideration. 



2. Kinetic description 

If a plasma contains magnetic irregularities that effectively 
scatter high energy particles in the vicinity of a shock front 
then these CR can be accelerated to the point where they have 
an important infiuence on the shock structure itself (the reader 
is referred to the reviews by Drury 1983, Blandford and Eichler 
1987, Berezhko and Krymsky 1988 and Jones and Ellison 1991 
for the details of shock acceleration not discussed here). With 
self consistent wave generation ignored, i.e. the spectrum of 
fiuctuations causing the scattering is assumed given, the ther- 
mal component is described by the three conservation equa- 
tions for an ideal fiuid, containing the added force exerted by 
the CRs. The transport equation for the isotropic part of the 
energetic particle phase space density, f{x,p,i), describes the 
CR interactions on a microscopic level. In one spatial dimen- 
sion these relations become 
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where p, U, Eg and Pg denote the fiuid's mass density, mass 
velocity, internal energy density and pressure with the time 
derivatives on the left hand side taken in a frame comoving 
with the fiuid. With a non-relativistic equation of state for the 
gas, Pg and Eg are related by 

Pg = (tg - 1)-Bg 

where 70 = 5/3. Equation (2.4) describes how the CR distri- 
bution function changes as a result of adiabatic compression 
and scattering off the magnetic turbulence. The CR pressure 
and energy density are given by appropriate moments of the 
calculated spectrum 
4ir 
~3 
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j P^i'f{x,p,t)dp 
Ec{x,t) = 4T j p^T{p)f{x,p,t)dp 



with T(p) the kinetic energy of a particle of momentum p. The 
local time dependent CR adiabatic index is then determined 
by 7c{x,t) = 1+Pc/Ec. 

This set of coupled partial differential equations can only 
be solved numerically. We adopt the method of Duffy (1992) 
which uses a finite difference scheme for equations (2.1)-(2.4). 
The numerical details are not repeated here. Physically we con- 
sider a piston driven shock with no initial CRs present. The 
process of injection whereby thermal particles are given an ini- 
tial boost to suprathermal energies, where their gyroradius is 
larger than the shock thickness, is not yet understood and is 
therefore modelled in the following, simple manner. A frac- 
tion, T], of the upstream gas particles incident on the shock per 
unit time are injected as CRs. These fresh CRs have a single 
momentum pi which is fixed at a value of 0.1 mc throughout 
this paper with m and c standing for the rest mass of an in- 
dividual (single species) CR particle and the speed of light re- 
spectively. The injection profile is therefore a delta function in 
phase space; at pi and the instantaneous position of the shock. 
Once the CRs appear they diffuse back and forward across the 
shock thus gaining energy by compression. The momentum 
dependence of k determines the range of diffusion lengthscales 
(i ~ k/U) and acceleration timescales (r ~ k/U ) from pi up 
to the highest energy particles. It would therefore be a good 
test of any model of the spectrum dependent closure parame- 
ters to choose a /c(p) that varies strongly with momentum since 
in this case both L and r will differ by several orders of magni- 
tude over the energy range of the CR spectrum. To this end we 
address that limit of CR transport theory (Skilling 1975) where 
strong wave excitation results in significant particle scattering 
after each gyration in the mean magnetic field. Thus the mean 
free path A is equal to the particle gyroradius; the Bohm limit. 
This gives k oa pv so that the diffusion coefficient scales like p 
(p) for non-relativistic (highly relativistic) particles. 

Before presenting a typical calculation it is worth noting 
briefiy how this system behaves in the, analytic, test particle 
limit; i.e. with the CR pressure gradient omitted from equa- 
tion (2.2). With the initial conditions described above, a pis- 
ton moving supersonically into an undisturbed and uniform 
medium, the gas dynamic solution is simply one of a shock 
wave that moves away from the piston separating uniform up- 
stream and downstream states. A power law CR spectrum 
f{p) = fo(plpi)~'' is produced at the shock from injection up 
to some time dependent cut-off pmax 



dt 
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with the acceleration timescale given by 

r = ? f ^ + 

" Ui - U2 \Ui ^ U2J 

where fiow velocities in the shock rest frame are denoted by U 
and the subscripts 1 and 2 refer to the upstream and down- 
stream regions respectively. The spectral index is related to the 
shock compression ratio, r = Ui /U2, by g = 3t/[t — 1) and the 
amplitude of the distribution function at pi is given by 
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with pi the upstream mass density. While the downstream CR 
distribution is uniform, and equal to that at the shock, the 
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Fig.2 Solution of the kinetic code for Mj, = 27.1 and j; = 0.002. 



upstream spectrum falls off exponentially with a lengthscale 
k/U. 

With the reaction included, the CR pressure gradient in 
the upstream region accelerates the plasma ahead of the shock 
to give a weaker subshock led by a broad precursor region. In 
figure 2 the first two panels show the velocity and pressure 
profiles for a piston Mach number of 27 and r] = 0.002. These 
plots are taken after the bulk fluid properties have reached a 
quasi-steady state where, although the CR spectrum at the 
high energy end still continues to evolve very slowly, the en- 
ergy densities do not change appreciably (Falle and Giddings 
1987, Duffy 1992). The evolution of the first, numerically cal- 
culated, closure parameter is shown in the third panel, with to 
the test particle acceleration timescale at injection. Initially fc 
(at the shock) is close to 5/3 since the only particles present 
are those with momenta close to pi. The CR adiabatic index 
decreases thereafter as the proportion of relativistic particles 
increases. Eventually, since the acceleration timescale goes at 
least linearly with momentum, it becomes harder to increase 
the upper cut-off and fc flattens off. 

The fourth plot shows the final CR spectrum at the shock. 
It is normalised throughout to fs{p) = fo(plpi)~''' which is the 
test particle spectrum that would be expected from subshock 
compression alone; fo is given by equation 2.8 applied to the 
subshock and = 3ts/{ts — 1) with Ts the subshock compres- 



sion ratio. Consider first those particles near injection. They 
have relatively small values for L and r. Consequently, as the 
initially strong shock evolves to give a weaker subshock plus a 
precursor, they are accelerated sufficiently quickly and over a 
small enough spatial scale to feel only the compression due to 
the evolving subshock. The spectrum near injection is therefore 
given, in amplitude and slope, by fs{p) which is evident from 
the normalisation used in the figure. Although perhaps not un- 
expected that test particle theory should apply to the lowest 
energy particles near the discontinuity it is by no means a triv- 
ial result and confirms the predictions of the double power law 
model of VDDM for the lowest energy particles. The spectrum 
is not so simple towards the higher end; here a single power law 
does not exist. The shock and precursor region evolve on the 
same timescale as that of the CR pressure which is determined 
by some intermediate momentum between pi and the upper 
cut-off. Since the highest energy particles are accelerated on a 
longer timescale than that of the fluid's evolution it follows that 
they cannot, unlike the CRs near pi, fully adjust to the instan- 
taneous fluid profile. It is therefore apparent that the shape of 
the spectrum at the highest energies is determined by the entire 
history of the subshock and precursor. These arguments can be 
tested by comparing the spectrum shape with that expected 
from the double power law argument qt = 3Tt/{Tt — 1) (dashed 
line). Although the spectrum flattens towards higher energies 
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the slope is nowhere quite as flat as might be expected from the 
total subshock-precursor compression. In addition the rollover 
of the spectrum towards the highest energies is not described 
by a cutoff with a qt spectrum (dot-dashed line). While this 
was expected it further weakens the quantitative applicability 
of a double power law approximation so that for strong shocks 
with a strongly momentum dependent k the double power law 
model will not suffice as a tool for calculating the closure pa- 
rameters in the fluid limit; the double power law spectrum will 
always overestimate the amount of highest energy CRs and will 
therefore give a lower value for fc ■ As an alternative it appears 
better to describe the overall spectrum by a power law at low 
momenta plus a bump at high momenta. We shall now con- 
struct a model for fc and 7c that is consistent with the generic 
features of such spectra. 
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3. Two fluid model and the closure problem 

The CRs can be quantitatively described as a fluid by taking 
the energy density moment of equation (2.4) to give 
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The adiabatic compression term now depends on the CR equa- 
tion of state through fc - The equations for Eg and Ec, (2.3) 
and (3.1), differ through the diffusion term for the CRs. The 
mean coefficient 7c is then related to /c(p) and the gradient of 
the spectrum by 



k{x) j p^T{p)^dp = j /CI 
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While we have gained, through the introduction of (3.1), a sim- 
pler description of the coupled system this has been done at 
the expense of introducing the two undetermined closure pa- 
rameters. With the numerical solution of the previous section 
as motivation the problem is now to construct an argument for 
the evolution of fc and 7c in terms of only those variables that 
enter the two fluid model. 

With a view to modelling integral moments of the CR dis- 
tribution, as opposed to predicting the exact distribution itself, 
our scheme essentially consists of approximating the power law 
and bump shape of a self consistent spectrum as shown in figure 
2 with a power law and spike. In other words the distribution 
function at the shock is written as 



/(p,<) = /o(^) +fiS{p-pt) 



The first term is the subshock power law distribution which 
is taken to extend from pi to the upper cut-off pmax(t) ob- 
tained by integrating equation (2.7) between the far upstream 
and downstream states over the fiuid's history and based on 
the total compression ratio qt. The effect of the enhanced dis- 
tribution at high energies is contained in the delta function 
taken at pb which is some constant fraction of the upper cut- 
off, a = Pb{t) / Piaa.xi{t) , the one free parameter in this approach. 
While equation 3.3 may appear rather primitive it contains all 
the predictable elements of the low energy part and makes only 
the simplest assumptions about the high energy spectrum. As 
a model for the spectrum of the highest energy particles equa- 
tion (3.3) is quite crude but we are more interested here in 
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predicting the closure parameters which are integral quantities 
as far as the momentum dependence is concerned. 

We can immediately write down zeroth order approxima- 
tions to the CR energy density and pressure as predicted from 
the test particle subshock spectrum 
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In a two fiuid calculation either equation (3.1) or one like it is 
solved to give Ec. The difference, Ei = Ec — Eo, is therefore 
the energy density of the extra component to the spectrum at 
high energies. By assumption in equation (3.3) this addition is 
entirely due to the delta function contribution near the upper 
cut-off so that the corresponding correction to the zeroth order 
pressure is 



Pi = (71 - 1)-Bi 



(3.6) 



where 71 = 1 + PbVb/3T[pb) is the adiabatic index of a delta 
function distribution at pb. An alternative but equivalent way 
of looking at this is that Ei determines the scale height fi by 



(3.3) fi 
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from which Pi can also be calculated. Either way we may now 
write down a model adiabatic index at the shock 
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Po+Pi 
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(3.8) 



that depends only on quantities which enter a two fiuid model 
and containing only one free parameter a of order unity. 

Turning to the mean diffusion coefficient it is not clear 
whether it is the most natural way to approximate 7c at the 
shock and then use this value at all positions in order to solve 
the system of equations (2.1)-(2.3) and (3.1). Since 7c does not 
connect local quantities (like fc with regard to Pc and Ec) 
but rather describes the diffusive character of the energetic 
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particle distribution, a presumably more appropriate way is to 
use a spatial average for 7c over the shock precursor in equation 
(3.1) to a first approximation, VDDM. To obtain such an aver- 
age we spatially integrate eq. (3.2) from the subshock at a; = 
into the upstream medium over a distance L where f{x,p,i) 
and Ec(x,i) go to zero 



r . -dEc ^ r, f 

I ax K— — = 4ir / ax I I 
Jo Jo J 

Defining the average < 7c > by 



= 4x / dx p'dpT{p)K^ (3.9) 
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and the precursor average < /c > by 

j dx j p^ dpT{p)K^ = j p^dpT{p) < K > f{x = Q,p,t). 

Combining the above we obtain the required precursor aver- 
aged mean diffusion coefficient as the energy weighted mean of 
< /c> 



< /c> = 



7 



4ir / p^dpT{p) <K> fs{p,t) 



4x / p^dpT{p)fs{p,t) 



7 



(3.10) 



For the model spectrum of equation (3.3) and the diffusion 
coefficient written as /c(p) = g{p)Ko , where /co is the value at 
injection, the second closure parameter then becomes 
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where dEo = 'lT^foP^T[p)[p/pi)~^' dp is the partial energy den- 
sity of the zeroth order spectrum. 

The coupled thermal and non-thermal system can now be 
solved numerically as a two fiuid system using equations (2.1) 
to (2.3) for the gas coupled to equation (3.1) for the CRs with 
the closure parameters derived from the model spectrum (3.3) 
as described above. There is one remaining problem of what 
exact value ought to be taken for a, the ratio of the Pb{t) to 
Pmax(t), or indeed even if it is sensible to tie the evolution of 
the bump to that of the upper cut-off by taking a to be con- 
stant. Although this can only be settled by a direct comparison 
between the results of a two fiuid model and kinetic solution 
(which is carried out below), the choice of a can be made less 
arbitrary by using the results of the full kinetic solution of the 
previous section. Here we can calculate the excess energy den- 
sity in the modified distribution over the test particle subshock 
distribution, AE = Ec — -Bo, and similarly AP = Pc — Po- 
These differences are due to the fiatter distribution towards the 
upper end of the spectrum in figure 2. We then ask at what 
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Fig. 5 Evolution of Pg and Pc for the fluid (solid) and kinetic (dashed) solutions with a range of Mach numbers and t} 

I < K{x,t) > 



0.002. 



Mj, 
17.81 
23.25 
27.10 
31.97 



7m /tC 

0.98 
0.98 
0.98 
0.97 

table 1 



2.4 
2.9 
3.7 
3.4 



momentum, p', we would have to place a delta function distrib- 
ution to obtain the correct ratio of AP/ AE. The amplitude of 
the delta function can also be obtained but is of little practical 
interest to us since it is not needed for the closure parameters. 
We then have 



p V 

2V) 



AP 
AE 



The evolution of ai[inetic(t) = p'(*)/Pmax(t) is plotted in figure 
3 for four different Mach numbers. The slow decline of akinetic, 
which in part can be attributed to the fact that pmax(t) as 
calculated by (2.7) is an increasing overestimate for CR spectra 
at modified shocks, suggests a value of a = 0.6 — 0.8 for the 
model spectrum (3.3). 

The numerical solution of the two fluid system is compared 
with that of the kinetic code in figure 4 for a piston Mach 
number of 23.25 and r] = 0.002 and for what is found to be the 
best value of a = 0.78. 



In each panel the solid and dashed lines show the evolution 
of quantities obtained from the fluid and kinetic codes respec- 
tively. The first two plots show the fall and rise of the thermal 
and non-thermal components of the total pressure at the gas 
subshock. The third panel compares the CR adiabatic expo- 
nent as given by equation (3.8) with that obtained from the 
CR spectrum at the shock with the two differing by no more 
than 2% throughout the evolution. The precursor averages of 
7c are shown in the fourth panel and differ by about a factor 3. 
What is most striking is that by replacing the local quantities 
fdxjt) and 7c(a;,t) with the spatially independent variables 
7m(t) and < /c(t) >m (eqns. 3.8 and 3.11), we have come close 
to reproducing the true, kinetic evolution of the bulk pressures 
with a two fluid model. This agreement is also found for a 
range of Mach numbers, at the same level of injection and with 
a = 0.78, as shown in figure 5. Of the two closure parameters, 
the time asymptotic value of 7m always lies within 3% of the ki- 
netic value at the shock which is a better agreement than that 
found for the precursor averaged diffusion coefficients which 
differ by a factor of between 2 and 4 (table 1). Given that we 
are comparing a two fluid model containing spatially constant 
closure parameters with a kinetic solution where fc and k are 
local quantities, this level of disagreement is perhaps not too 
surprising. Figures 6a and 6b show how, for Mp = 23.25, small 
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Fig. 6a Pg from the kinetic (solid) and fluid codes for a = 0.75 
(dot-dashed), 0.78 (dashed), 0.8 (dotted). 

changes in the parameter a alter the evolution of the two fluid 
pressures. 

4. Conclusions 

We have used a model of CR spectra in time dependent modi- 
fied shocks to predict values for the closure parameters of CR 
hydrodynamics at shock fronts. The lowest energy particles 
have a spectrum determined by the structure of, and injection 
at, the gas subshock. This accounts for almost all of the CR 
energy density at the subshock, the excess is due to more ener- 
getic particles that feel the extra compression in the precursor 
and have a flatter spectrum. We have looked at intermediate 
timescales where the most energetic CRs, whose acceleration 
timescale is so large that the slope of their spectrum cannot 
adjust to that given by the double power law model. For our 
dynamical, trans-relativistic, range (chosen as an intermediate 
case where 7c is not close to either of its limits) it is more 
accurate to attribute the excess energy to an individual parti- 
cle momentum which evolves with the upper cut-off. The full 
hydrodynamic description of the coupled interaction between 
the thermal gas and the non-thermal CRs then agrees quite 
closely with that of a full kinetic treatment. As mentioned in 
the introduction there is evidence from a recent kinetic solu- 
tion (Berezhko et al. 1993) that, for piston driven shock so- 
lutions, with an increased dynamical range the CR spectrum 
does eventually relax to the original double power law model.* 
For timescales much longer than those considered here, where 
7c eventually evolves towards 4/3 in the planar case, such a 
picture could then be used to calculate the closure parameters. 
Acknowledgements We would like to thank J. Kirk and 
E. Berezhko for discussions. 
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* Ultimately the rate of particle acceleration at the shock 
front of a SNR will decay with decreasing Mach number to 
give low and high energy parts of the spectrum from particles 
accelerated during this final phase and older CRs that have 
been accelerated earlier. 
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Fig. 6b Pc from the kinetic (solid) and fluid codes for a = 0.75 
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